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We study simple lattice systems to demonstrate the influence of interpenetrating bond networks 
on phase behavior. We promote interpenetration by using a Hamiltonian with a weakly repulsive 
interaction with nearest neighbors and an attractive interaction with second-nearest neighbors. In 
this way, bond networks will form between second-nearest neighbors, allowing for two (locally) dis- 
tinct networks to form. We obtain the phase behavior from analytic solution in the mean-field 
approximation and exact solution on the Bethe lattice. We compare these results with exact numer- 
ical results for the phase behavior from grand canonical Monte Carlo simulations on square, cubic, 
and tetrahedral lattices. All results show that these simple systems exhibit rich phase diagrams 
with two fluid-fluid critical points and three thermodynamically distinct phases. We also consider 
including third-nearest-neighbor interactions, which give rise to a phase diagram with four critical 
points and five thermodynamically distinct phases. Thus the interpenetration mechanism provides 
a simple route to generate multiple liquid phases in single-component systems, such as hypothesized 
in water and observed in several model and experimental systems. Additionally, interpenetration 
of many such networks appears plausible in a recently considered material made from nanoparticles 
functionalized by single strands of DNA. 

PACS numbers: 64.70.Ja, 64.60. De, 68.35.Rh 



I. INTRODUCTION 

In the last 15 years, liquid-liquid phase transitions in 
one-component systems has been an area of vigorous re- 
search [TJ [21 Liquid-liquid phase transitions have 
been found experimentally in phosphorus [3], and are 
also suspected for many tetrahedrally coordinated flu- 
ids including water [5J |BJ E EE US], carbon [TT], sil- 
ica [T21 H3] , and silicon [T4] . A variety of different ap- 
proaches have been used to understand the emergence 
of a second liquid state, most of which rely on a com- 
petition between non-directional van der Waals inter- 
actions and directional bonding interactions that favor 
more open states [L5l [TBI HZ1 EH EH] ■ Some success has 
also been found for models with purely symmetric inter- 
actions [2Q1 [2D 1221 H3] • Recently, a model for nanoparti- 
cles functionalized by single strands of DNA [2H [23 [2S] 
whose sequence promotes bonding between the nanopar- 
ticle units showed that polyamorphic behavior can also 
arise due solely to directional interactions that result in 
open networks which interpenetrate, giving up to three 
critical points and four amorphous phases |27j . The ques- 
tion that remains is, can network interpenetration be a 
general mechanism to liquid-liquid phase transitions? 

Lattice models have long served to provide analytic 
insight into complex physical problems [28 1 221 ISO] ■ The 
most famous example is the Ising model, which helped 
to understand the origin of spontaneous magnetization 
in magnetic systems. The Ising model can also be re- 
cast as a lattice gas model, which equivalently provides 
an understanding of condensation of the liquid from the 
gas. The Ising model (or lattice gas) exhibits such com- 
plexity while only using very simple first neighbor inter- 
actions. Making these interactions or lattice occupancy 



state more complex - for example, the Potts model, the 
spherical model, etc - can yield even richer behavior |28j . 

Indeed, lattice model approaches based on a lattice 
gas with additional orientation dependent bonding inter- 
actions [T51 [T71 HS1 EH] have reproduced phase behavior 
with two critical points. Motivated by the finding that 
interpenetration might serve as a mechanism for gener- 
ating liquid-liquid transitions, we explore lattice mod- 
els that incorporate interpenetration to see how phase 
behavior is affected. In our models we promote inter- 
penetration of networks by using an attractive interac- 
tion with second-nearest neighbors (2NN) and a weakly 
repulsive interaction with nearest neighbors (NN). The 
properties of Ising and lattice-gas models with NN and 
2NN interactions were first examined on a two dimen- 
sional lattice, primarily for the case of purely antifcr- 
romagnetic interactions (i.e. both NN and 2NN repul- 
sions) [3U 1321 (Ml El [Ml M M\- Our case (repulsive 
NN and attractive 2NN interactions) has been exten- 
sively studied for the 2D triangular lattice [3SJ EH1 HO] , 
where the results have been applied to gas adsorption on 
surfaces. Considerably less attention has been given to 
3D lattices; most work has focused on the zero-field Ising 
formulation of the model on the cubic lattice [HJ H2J E3] > 
where both the phase behavior and critical properties 
have been examined. The phase behavior for the lattice- 
gas formulation has also been examined on the hexagonal 
close-packed lattice with application to binary metal al- 
loys [44]. Here we examine the lattice-gas formulation 
of the model on 3D lattices including the cubic lattice 
and the diamond lattice. In the Ising formulation, this 
is equivalent to including a non-zero external field. For 
the completeness of our work, we also consider the be- 
havior on the 2D square lattice. For all three lattices 
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we carry out both the mean field solution and exact nu- 
merical Monte Carlo (MC) solution. For the symmetric 
lattice model, we also provide an exact solution on the 
Bethe lattice. 

Our findings indicate that on all lattices there are two 
critical points and three distinct phases: (i) unassociated 
molecules of gas, (ii) liquid I, a single network of alter- 
natingly filled sites, and (iii) liquid II, a double interpen- 
etrating network with all sites occupied. For all lattice, 
the two critical temperatures are the same due to the 
symmetry between occupied and empty sites. To con- 
nect more closely with systems of experimental interest, 
we also consider a slightly more complicated model that 
includes a three-body term that accounts for increased 
repulsion in crowded states. This model breaks the sym- 
metry between occupied sites. For this model, we find 
that the high density critical point occurs at a tempera- 
ture lower than the low density critical point, as observed 
in most polyamorphic systems. 

The Bethe lattice has a distinct geometry from other 
lattices we study, and it exhibits somewhat different 
phase behavior. While there are still two critical points 
and three distinct phases, there is at least one case 
in which the coexistence lines merge for temperatures 
slightly lower than T c and separate again at even lower 
temperatures. In addition, a density anomaly occurs, as 
does in water |5]. 

Since it has been argued that the observation of even 
more amorphous phases can result from very open net- 
work structures [27], we also consider a modification of 
the previous model. Specifically, we consider a model 
including third-nearest-neighbor (3NN) attraction, and 
weak first and second neighbor repulsion [33 . Results on 
the square lattice show that there can be up to 4 critical 
points for this system. 

The paper is organized as follows: Sec. [TTJdescribes the 
three models: (i) the symmetric model, (ii) the asymmet- 



ric model, and (iii) the 3NN interaction model. Sec. Ill 



gives the derivation of the mean-field approximation and 
Bethe lattice solution, and describes the methods used 
in the MC simulations. Sec. |IV| presents the results and 
discussion. We conclude briefly in Sec. |V] 



the total number of sites taking the volume of a single 
site v = 1. We define the ratio of interaction strength 
R=e 1 /e 2 . 

To promote interpenetration, we choose ei < and 
(-2 > 0. In this way, bond networks will form between 
second-nearest neighbors, allowing for two (locally) dis- 
tinct networks to form. In nature, this phenomenon oc- 
curs in ices VI, VII and VIII, where two interpenetrating 
structures form. Multiple interpenetrating networks also 
occur in model of DNA-functionalized nanoparticles, and 
gives rise to polyamorphic phase behavior [27] . The aim 
of our lattice model is to retain this interpenetration fea- 
ture while eliminating other complexities. We expect this 
simple lattice model to exhibit multiple critical points 
and to serve as a demonstration of the influence of inter- 
penetrating bond networks on phase behavior. 

This lattice model (like the Ising model) has an intrin- 
sic symmetry between occupied states and unoccupied 
states; thus we call this the symmetric model. 



B. Asymmetric Lattice Model 

The intrinsic symmetry between occupied and unoc- 
cupied states is usually not found in nature. To break 
the particle/hole symmetry, we introduce a three-body 
term in the NN interaction. Such three-body interac- 
tions for the antifcrromagnetic case have been studied 
elsewhere [45] [46] . This term accounts for the increased 
repulsion as crowding occurs. One can think of this three- 
body interaction as scaling the NN interaction strength 
by the ratio of number of occupied sites around the pair 
to the number of sites around the pair. This leads to the 
Hamiltonian 



NN <ijk> 



IJ_k_ 

C 



62 Y nin j' 
2NN 



(2) 



where J2<ijk> sums over a U sites k that are first neighbor 
of either site i or site j; c is the total number of such 
sites for any pair On the 2D square lattice and 3D 

diamond lattice, c = 6; on a cubic lattice c = 10. 



II. MODELS 
A. Symmetric Lattice Model 

The first and simplest model we consider is a second 
neighbor lattice gas with Hamiltonian 

H = -ei} njnj - e 2 )] rijUj , (1) 

NN 2NN 

where J2nn indicates a sum over all NN pairs and J^2NN 
indicates a sum over all 2NN pairs. At each site i, the 
occupancy is 1 when the site is occupied and when 
the site is unoccupied. The volume V of this system is 



C. Third-neighbor Interaction Model 

We shall see that the lattice model with only first and 
second neighbor interactions cannot account for more 
than two critical points. Hence, as an extension, we con- 
sider a model including up to third neighbor interactions. 

T-i = ~ei ^ iiiUj - e 2 Y n i n j ~ e 3 Y n%n i ( 3 ) 

AW 2NN 3NN 

As in previous models, we choose t\ < 0, £2 < 0, and 
£3 > to promote interpenetration. We define two ra- 
tios of interaction strength: R\ = ei/e^ and R2 = £2/63. 
This model is more representative of systems in which the 
bonding is very large in comparison to the core repulsion. 
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FIG. 1: Illustration of sub-lattice division for the MF ap- 
proximation, (a) 2D square lattice, (b) cubic lattice, (c) 
diamond lattice, (d) 2d square lattice with 3NN interaction. 
For (a)-(c) the division is made such that there is no NN or 
2NN interaction between sites on the same sub-lattice. For 
(d) the division is made such that there is no NN, 2NN, or 
3NN interaction between sites on the same sub-lattice. 

This third neighbor model for the purely antiferromag- 
netic case (ei < 0, £2 < 0, £3 < 0) has also previously 
been examined [33 . 

III. METHODS 
A. Mean Field Approximation 

1. Symmetric Case 

To carry out the mean- field solution, we follow the pro- 
cedure described by Binder and Landau |32j . Specifically, 
we first divide the lattice into sub-lattices such that sites 



TABLE I: List of relevant values for the MF approximation. 
f(a,/3) and g(a, j3, 7) are used for the asymmetric model. 
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on the same sub-lattice do not interact with each other. 
The way we make such division for the three lattices is 
illustrated in Fig. [I] We denote the number of such sub- 
lattices v. We define the interaction energy between a 
site on sub-lattice a and a site on sub-lattice (3 as e a p. 
We denote the number of such (a, (3) pairs for a site on 
either sub-lattice as 70/3. If the a and (3 sub-lattices are 
nearest neighbors, e a p — t\ and 7 Q( g = 71. If a and (3 
sub-lattices are second-nearest neighbors, e a p — 62 and 
lafi = 72- Otherwise there is no interaction so e a p — 
and 7 Q( g = 0. The values for v, 71, and 72 for the three 
lattices we study are listed in Table [I] 

In mean-field approximation, the occupancy of a neigh- 
boring site is approximated as the mean density of the 
sub-lattice to which it belongs. There are V/v such sites 
on each sub-lattice, thus the mean- field Hamiltonian is 

H MF = — — ^2 7a/3^a[3PaPl3 = ~y ^ la^aP^aN 'p, 
a,/3 a,/3 

where p a = {N a v)/V, and N a — ^2 i£a rii is the number 
of occupied sites on sub-lattice a. To avoid double count- 
ing, the sum Y, a ,p is specifically Ea=iEL a+ r Tnc 
mean-field Hamiltonian depends on the configuration via 
Ni : N 2 ,...N v . We label this set {N a } for brevity. 

We fix {N a }, V, T, and evaluate the canonical parti- 
tion function 



Z({N a },V,T) = 



n 



V/i 



a,0 



(5) 



r 



Using the thermodynamic relations and Stirling's approximation, we find 

Pa 



li a = kT\TL—^ ^2j a f3e a ppp, (7) 

1 Pa 13 



\0N a J {N ^ a}VT \ dV J {N a },T P = ^ln(l - p a )- -^j a pe a pp a pp. (8) 

(6) V a V a,0 
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Ref. 32J derives this result using the grand canonical particularly useful when we introduce three-body terms, 
partition function. Using the canonical approach will be Rearrangement of Eq. Q leads to 




where a goes from 1 to v. Eq. ([9]) represents an system of v implicit equations for sub-lattice densities {p a }. Given 
p and T, we obtain the set {p a } by solving this system of equations numerically using the iterative scheme described 
in Ref. [32]. Specifically we use 



P^ 



(cos 2 0)pi"- 1) +(sin 2 



jexp 



y^^laptappp - Pa /kT 



1 + exp 



> (io) 



where <fi is an arbitrary tuning parameter. When = 0, 
the series stays at p a °^- When <f> = ir/2, p a is a direct 
iteration of Eq. ^ . One needs to choose appropriate 
so that the series converges. As discussed in Ref. [52] , 
this iterative scheme prevents us from finding unstable 
solutions which are not physical. We choose the initial 
densities {p a } and iterate using Eq. (10 1 until — 

pi™ -1 '' I < S for all a. In our calculations the tolerance 
5 = 10~ 9 , and most of the time we choose = 7r/4. 

Eq. ^ can have multiple solutions. To find all possi- 
ble stable or metastable solutions, we repeat the iteration 
procedure with different sets of initial densities. In the 
cases when multiple solutions are found, we calculate the 
pressure using Eq. (|8|, and take the solution with the 
highest pressure (i.e. lowest grand potential free energy 
f2), which is the most thermodynamically stable state. 
We can distinguish first order transitions by a disconti- 
nuity in density (i.e. and second order transitions 

by a discontinuity in the slope of density (i.e. §~t). 



write out the NN and 2NN terms separately, unlike the 
symmetric case. The mean-field Hamiltonian is 



ftMF _ __ 



V 



7i£i 
c 



22,p<xPf3 ^2 ps 

NN <af3S> 



72^2 ^2 paf) P 
2NN 

(11) 

where ^2 NN sums over all distinct pairs of sub-lattices 
(a, (3) that are first neighbors, ^2 2 nn sums over all dis- 
tinct pairs of sub-lattices (a, 0) that are second neigh- 
bors, and X)< Q /3i5> sums over all c sub-lattices 5 that are 
a first neighbor of sub-lattice a or (3. It is possible that 
a particular sub-lattice may appear twice as a neighbor 
to lattice a and /3, so some sub-lattices must be counted 
twice in this summation. 



2. Asymmetric Case 



For the asymmetric model, the NN and 2NN interac- 
tions have a different form from each other, and we must 



Now rewrite the Hamiltonian as a function of {N a }, 



(12) 



NN <aP&> 

and write out the canonical partition function, fixing {N a }, V, and T, 



2NN 



Z({N a },V,T) 



n 



V/v 
N a 



exp 



NN 



<a/38> 



2NN 



/kT 



(13) 
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To calculate /i, we take partial derivative of the summation J2nn with respect to N a . When doing so, we need to 
consider that N a can appear in the NN sum, as well as the three-body sum, so there will be 3 terms associated with 
the derivative of the near neighbor interactions. Accordingly, we obtain 



[i a = fcTln 



1 



7i£i 
c 



J2 

a.NN 



^2 + Pa f (a, (3) 



r 



E 

NN 



a,2NN 



P/3, 



(14) 



9E ir;; >w , and <?(<*, /?,<*) = 

-g- . The summation nn sums over an P 

that are first neighbors of a, J] jvjv 



where /(a, (3) 



sums over all 

pairs (P, d) where (3 and 5 are first neighbors of each 
other, but (3 ^ a and 5 ^ a, and J2 a 2NN sums over all 
(3 that are second neighbor of a. 

f(a,/3) quantifies the repulsion on a due to 3-body 
terms containing only a and (3, and g(a, (3, 8) quanti- 
fies the repulsion on a due to 3-body terms containing 



a, (3, and 5. For the 2D square lattice, f(a,(3) = 1 
and g(a,(3,5) = 2. For the cubic lattice, f(a,(3) = 1, 
and g(a,f3,5) can be or 2 depending on weather a is 
a neighbor to the pair (f3, 8). For the diamond lattice, 
f(a,(3) = and g(a,(3,5) = 1. These values are also 
listed in Table U 

To calculate P, we take partial derivative of In Z with 
respect to V. Here we get an extra factor of 2 for the 
NN interaction term because of the power (v/V) 2 , which 
arises from its three-body nature. We obtain 



kT 



2 7i£i 

2^ Pap P 



cv 



PS 



NN 



<a/3S> 



72^2 \ ^ 
— 2^ P«PP 
2NN 



(15) 



Similar to the symmetric case, Eq. ( 14 1 is a system of v 



implicit equations for {p a }- We rearrange it into the form 
of Eq. (|9| to carry out the iteration scheme of Eq. (10 1. 



We repeat this iteration procedure with different sets of 
initial densities, and take the solution with the highest 
pressure (lowest f2). 



B. Solution on Bethe Lattice 

While exact solution on most lattices is cither difficult 
or impossible, the recursive nature of the Bethe lattice 
(which we define explicitly below) sometimes makes so- 
lution possible. The Bethe lattice usually offers a better 
approximation to a regular lattice with the same coordi- 
nation than the conventional mean-field approximation, 
since it preserves correlations between local sites. 

A Cayley tree, as illustrated in Fig. [2] is constructed 
as follows: 

(i) We start from a center site and connect 7 sites to 
it. We say the center is generation M, and the 7 
nearest-neighbor sites are generation M — 1 . 

(ii) We connect 7— 1 sites to each site in generation M— 
1. Repeat this procedure until we reach generation 
1. 

For most lattices, the ratio of surface sites to interior 
sites vanishes in the thermodynamic limit M — > 00. How- 




FIG. 2: Illustration of a Cayley tree with 7 = 3, 
Generation number is labeled on the graph. 



M 



5. 



ever, on a Cayley tree, this ratio does not vanish [28] . 
The standard way to overcome this problem is to con- 
sider only sites that are infinitely far from the boundary 
in the thermodynamic limit. These central sites will not 
be affected by the surface, and they constitute the Bethe 
lattice [2S]. 

We now describe the solution to the symmetric lattice 
model on the Bethe lattice. We use um to denote the 
occupancy for the center site and {n m } to denote the 
set of occupancy for sites in generation m. We use 3m 
to denote the grand partition function of the whole tree 
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7m, and use j m to denote the partial grand partition appear explicitly, so 3 m is a function of n m and n m -i- 
function of a branch B m whose top site is in generation Specifically, 
to. When evaluating j m , the occupancies n m and n m -\ 



£ ex P 

{«*} 



13 if! ^2 n i 



ei £ r^n,- + e 2 £ n,^ 

NNeB m 2NNeB m , 



(16) 



where /3 = 1/kT, not to be confused with the lattice enu- 
merating variable in Sec. Ill A X){n } sums over all pos- 
sible configurations in branch £>,,„ , fixing n m and n m _!. 
SiGB sums over all sites in B m except for the top site in 
generation to. i m does not account for 2NN interaction 
between different branches of the same generation to. 

The whole tree Tm is composed of 7 branches Bm, so 
the total grand partition function 3 m can be written as 



Eq. (17 1. In our solution, we explicitly account for 2NN 



{iiu-i} 



'M-lJ ) 



(17) 

where ^ sums over all pairs (i,j) of sites in generation 
M — 1. The term e^ M ™ M takes into account the occu- 
pation of the center site, and the term e' 362 ^ ni ™ J takes 
into account the interaction between different branches 
due to the 2NN interaction between the 7 sites in gener- 
ation M—l. Fig.j3ja) gives a graphical representation of 



interaction between branches of the same generation us- 
ing the same approach as Ref. |47j . Such branch inter- 
actions do not appear in nearest neighbor only models, 
and many longer ranged models in the past have ignored 
branch interactions. 

When n of the 7 sites in generation M—l are occu- 
pied, the interaction energy of those sites across branches 
of the same generation is — £2(2)- There are Q) such con- 
figurations. Thus we can rewrite 3 m as 



3m - £ e^™« £ (A eM2) 3M (n M , 1) iV l (n M , 0) . 

(18) 

A branch B m is composed of 7 — 1 sub- branches B m -i, 
so the partial grand partition function j m (n m , n m _i) can 
be written as 



7-1 



1 "m-2/ ! 



(19) 



where sums over all sites i in generation m — 2, and 
^2 sums over all pairs (i, j) of sites in generation to — 
2. The term e P^ n ™-± takes into account the occupation 
of the site in generation to — 1, the term e P^ n mn m ^i 
takes into account the NN interaction between the site 
in generation m and the site in generations to — 1, the 

term e' 9e2 ™ m ^- ni takes into account the 2NN interaction 



between the site in generation to and sites in generations 

and to — 2, and the term e 13 ^^ nirlj takes into account 
the 2NN interaction across sites in generation to — 2 of 
separate branches. Fig. [3^b) gives an easier to follow 
graphical representation of Eq. (19 1. 

We can rewrite ]> m { n m, n m-i) using the same ap- 
proach used for Eq. (Il8| as 



1 E ( 7 n ^e^^^in^^f^rin^O). (20) 

n=0 ^ ' 



This is the recursion relation for our model on the Bethe lattice, which facilitates analytic solution. We consider 



iM 



FIG. 3: Graphical representation of the separation of the Bethe lattice grand canonical partition function into different branches. 
We represent NN interaction with solid lines, 2NN interaction with dotted lines, and indicate the contribution from the chemical 
potential with filled circles. Empty circles represent sites where the chemical potential contribution is not counted. The left 
and right panels of (a) corresponds to the left and right hand sides of Eq. (17 1. It shows how the total grand partition function 
5m can be divided into (i) 7 partial grand partition functions jm, (ii) 2NN interaction between the branches (represented as 
a triangle), and (iii) occupation of the center site (represented as a filled circle in the triangle). The left and right panels of 
(b) corresponds to the left and right hand sides of Eq. ( 19 1. It shows how the partial grand partition function j m can be dived 
into (i) 7 — 1 partial grand partition functions of shorter branch $ m -i, (ii) 2NN interaction between generation m and m — 2 
(represented as a triangle), (iii) occupation of the site in generation m — 1 (represented as a filled circle in the triangle), and 
(iv) NN interaction between the site in generation m and the site in generation m — 1 (represented as a solid line). 



the possible values of n m and n m —i of Eq. (20 1, which 
can be simplified if we define the following 

_ 3m (1, 1) , _ 3m (1,0) _ 3m (0, 1) 

am "3m(0,0) Om - 3 m(0,0) Cm "3m(0,0)- 

_ (21) 
We evaluate Eq. (21) using the four cases for n m and 
n m -i in Eq. (20 I. In the limit M — > oo, we obtain 



(22) 



TP 



-l) e M3)cn 
i / 



where a — limM— >oo am , b 
lim M ->oo cm- 



lim 



M- 



(24) 



and c 



From Eq. ( 18 1 we can directly derive the expectation 



(23) value of the occupancy of the center site 



< n M >= 
and in the limit of M going to infinity 



E„ M ef*»" El=o (iy e2 ®& (nu, 1) 3X7" (n M , 0) ' 



p = lim < tim >= t^v ~ 



Af^oo 



£ til U 

The free energy is calculated using the method Gujrati proposed 

„ . r t-1 



— = - lim 

kT 2 m^oo 



In 3m - Yl ln 3A/-i 



fc=i 



(25) 



(26) 



(27) 



Note that this expression was derived by Gujrati in the specific case of nearest-neighbor interactions only. With 



8 



the graphical representations in Fig. [3] we can show that 
this expression is also valid for models with 2NN interac- 
tions (if interactions across separate branches are not ig- 
nored) . Specifically, we can divide the whole tree 7~m into 
7 branches Bm, and further divide them into 7 x (7 — 1) 
shorter branches Bm—i- In the graphical representation 
we get 7+I triangles and 7 lines in this process. We 
can divide the 7 — 1 shorter trees Tm-i into 7 x (7 — 1) 
shorter branches Bm-i- In the graphical representation 
we get 7 — 1 triangles and 7 — 1 lines in this process . 
The difference is 2 triangles and 1 line, which represent 
the occupations of 2 lattice sites, 7 x (7 — 1) 2NN inter- 



actions, and one NN interaction. This is the free energy 
of two sites. Thus we can divide it by two to obtain the 
free energy per site [i.e pressure P). 



Eq. (27 1 can be simplified by substituting 3 m and 
3 m- 1 using Eq. ( 18 1, and factoring out powers of Im (0 , 0) 
and 3m-i(0, 0) to express the rest with the ratios 6m, 
cm, «m-i, bu-i, cm-i- Further substituting 3m(0,Q) 
using Eq. (p0|), we can cancel out the powers of 3m-i(0, 0) 



and express the right hand side of Eq. (27) entirely with 
the ratios a M , b M , c M and a M -i, b M -i, c M -i- 
we take the limit M — > 00, we obtain 



When 




We have now derived everything we need to evaluate 
the phase behavior. To evaluate p and P for a given p and 
T, we first obtain a, 6, c by solving the system of implicit 
equations (??) using the iterative scheme described in 
Sec. |III A| and use these values to calculate density p 
with Eq. ( 26 1 and pressure P with Eq. ( 28 ) . Similarly 



we repeat the iteration procedure with different initial 
a, 6, and c to search for possible stable and metastable 
states. When multiple solutions are found, we take the 
solution with the highest pressure (lowest free energy). 



togram re-weighting to make minor adjustments in T and 
p so that the distribution P(M) has same height for the 
two peaks corresponding to the two phases, and that the 
integral under each of the peaks, that is the probability 
for each phase, is equal. 

Data used to build the histograms are obtained from 
pairs of (N, E) values taken from 8 independent simula- 
tions, each running 5 x 10 6 MC steps per site. The (TV, E) 
data are taken every 10 MC steps for the entire lattice 
after an equilibrium of 5 x 10 5 MC steps per site. 



C. Monte Carlo Simulations 

To evaluate the exact phase diagrams, we perform 
Monte Carlo (MC) simulations in the grand canonical 
ensemble (fixed p, V, and T) - or GCMC for short gS]. 
For the 2D square lattice, our system size is 40 x 40 (1600 
sites). For the cubic and diamond lattices, our system 
size isl2xl2xl2 (1728 sites); for the diamond lattice 
this corresponds to 6 x 6 x 6 unit cells. Periodic boundary 
conditions arc implemented in all cases. 

We first locate the rough location of critical points by 
the appearance of a bimodal density distribution. To 
obtain accurate T c , p c , and p c , we use the fact that 
the model is expected to be in the Ising universality 
class Specifically, we study the order parameter 

AI = p — su and use the histogram re-weighting tech- 
nique [50J to make minor adjustments in T and p so that 
the order parameter distribution P(M) approaches that 
of the Ising universality class. A detailed description of 
this procedure is given by Refs. |50j . 

To evaluate the phase boundaries in the subcritical 
region, we perform a series of GCMC simulations with 
multi-canonical biased sampling [49 , 50 to allow us to 
efficiently sample both phases in a single simulation. To 
estimate the coexistence densities, we again apply his- 



IV. RESULTS AND DISCUSSION 
A. Symmetric Lattice Model 

To illustrate how the long-ranged attractions of the 
model Hamiltonian Eq. ([IJ affect phase behavior, we fo- 
cus on the phase diagram for the case of R = —1/2, i.e. 
second neighbor attractions twice that of first neighbor 
repulsion. We show the resulting p-T phase diagrams for 
the mean-field approximation and for MC simulations on 
square, cubic, and diamond lattices in Fig. |4ja,b,c). On 
all three lattices we find three thermodynamically dis- 
tinct phases analogous to: (i) unassociated molecules of 
gas, (ii) liquid I, a single network of alternatingly filled 
sites, and (iii) liquid II, a double interpenetrating net- 
work with all sites occupied. Hence, our system demon- 
strates that long-ranged attraction and short-ranged re- 
pulsion alone can generate multiple high-density phases 
via interpenetration. This is precisely the mechanism 
proposed for DNA functionalized nanoparticles [2"T] . 

By construction of the model, the p-T phase diagrams 
are symmetric about p=0.5. The mean-field approxima- 
tion correctly predicts the number of transitions and the 
qualitative shape of the first-order transition boundaries. 
As expected, the mean-field approximation overestimates 



9 





(e) 2 



tt 


— o— .... 


■ u 




Square 




0.5 1 1.5 




& 


0" 


Cubic 








0.2 






FIG. 4: Phase diagrams from the mean-field approximation 
and Monte Carlo simulations of the symmetric lattice model 
with R — —1/2. (a,d) square lattice. (b,e) cubic lattice. (c,f) 
diamond lattice. 



the terminal temperature of the first order transition 
and incorrectly predicts a second-order phase transition. 
While the pressure is not readily available from the MC 
simulations, Eq. ^ allows us to readily calculate the T-P 
phase diagram for each lattice (Fig.[3Jd,e,f)). In all cases, 
the slope of the coexistence lines in the T-P plane are 
positive, unlike the slope of the coexistence lines in wa- 
ter. For the coexistence lines, the slope (f^) = (the 
Clausius-Clapeyron relation) [29"] ; thus, interpenetration 
alone does not require an anomalous relation between the 
difference in entropy As and difference in volume Av of 
the two phases. 

Since the Bethe lattice geometry is distinct (i.e. no 
closed loops) from those studied for the MF and MC so- 
lutions, we present the results separately. To mimic the 
behavior of the diamond lattice studied in Fig. [4^c,f ) and 
known to occur in many polyamorphic fluids, we examine 
the case 7 = 4 (Fig. 5|. Moreover, both the Bethe lattice 
with 7 = 4 and the diamond lattice have 12 second neigh- 
bors, so the results can be expected to be similar. To our 
surprise, the phase behavior for R = —1/2 is qualitatively 
different from that of the diamond lattice. Specifically, 
the coexistence lines merge for 1.754 < kT j 62 < 1.886, 
resulting in two triple points in the T-P phase diagram. 
In other words, there are only liquid and gas states for 
1.754 < kT/e2 < 1.886, but just below and above there 



FIG. 5: Phase diagrams of the symmetric model on Bethe 
lattice, calculated from exact solution. (a,c) R = —1/2. 
(b,d) R = —3/4. The shaded regions indicate the location 
of anomalous p dependence, i.e. ap < 0. 
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FIG. 6: Isochores of P at constant fj, for the Bethe lattice 
R = —3/4 case. Density anomaly exists for regions where the 
slope of the isochores are negative; such regions are shaded. 
The slope of the isochores between the two boundaries all go 
to zero at kT/e2 = 1.194, as indicated by the dotted line. 



are three distinct states. We also study R = —3/4; in 
this case, there are two distinct phase transitions that do 
not merge. 

The presence of a negatively sloped coexistence line in 
the T-P plane indicates an anomalous ratio ^ < 0. This 
suggest there may be anomalous density dependence, i.e. 

that the isobaric expansivity ap = — ^ (^§j^j < 0. 

Since our solution does not easily allow us to hold pres- 
sure fixed, we instead check the thermal pressure coef- 
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FIG. 7: Phase diagrams from the mean-field approximation 
and Monte Carlo simulations of the asymmetric lattice model 
with R — —1/2. (a,d) square lattice. (b,e) cubic lattice. (c,f) 
diamond lattice. 

ficient "fy — {^f)y i since 7y = where kt is the 
isothermal compressibility. The sign of ap is determined 
by the sign of jy since > 0. We find that for the 
R = —1/2 case, the density is anomalous throughout the 
intermediate phase. Fig. [6] shows that for the R = —3/4 
case, the density is anomalous for kT ' je^ < 1-194 in the 
intermediate phase. In the phase diagrams in Fig. [5] 
density anomaly regions are indicated by the shading. 



B. Asymmetric and Third- neighbor Lattice Models 

We next examine the changes in the phase diagrams of 
the symmetric model (except for the Bethc lattice) when 
we include a 3-body interaction with first neighbor sites. 
For the same value R = —1/2, Fig. [7] shows that the 
three-body term breaks the symmetry in p-T phase dia- 
grams. While the width of the transitions is unchanged, 
the high density critical point has a lower critical temper- 
ature than the low density critical point. The depression 
of the second critical temperature is consistent with the 
observation in polyamorphous systems that the high den- 
sity critical point occurs at lower T than the liquid-gas 
critical point. Aside from the difference in T c , the phase 
diagrams in both T-P and p-T are qualitatively compa- 
rable. 




FIG. 8: Phase diagrams of the 3NN interaction model from 
the mean-field approximation and Monte Carlo simulations, 
with Ri = i?2 = —1/2- Done on a 2D square lattice. 

The third-neighbor lattice model has a longer range 
for bonding, which opens more nearby sites that can be 
occupied by molecules in separate interpenetrating sub- 
lattices. Our MF and MC results for the phase behavior 
of the case R\ = R2 = —1/2 (Fig. [8| confirm that the 
additional open sites allow for an even greater number 
of coexisting phases; specifically, we find up to five co- 
existing phases - more than previously observed in any 
molecular systems. These five phases correspond to (i) 
gas, empty lattice, (ii) a single network of one sublat- 
tice, (iii) two networks with two sublattices occupied, (iv) 
three networks with three sublattices occupied, and (v) 
four networks with the whole lattice occupied. Although 
the coexistence line of the highest density transition in 
the T-P plane indicates an anomalous ratio ^ < 0, by 
examining jy, we find there is no density anomaly. 

The middle phase consists of sites occupied alternately 
in parallel (superantifcrromagnctic) instead of sites oc- 
cupied alternately in diagonal (antifcrromagnetic). This 
can be explained by the ground state free energy of these 
two configurations. For an occupied site on a fully super- 
antiferromagnetic configuration, half of the NN sites are 
occupied and none of the 2NN sites are occupied. For an 
occupied site on a fully antiferromagnetic configuration, 
none of the NN sites are occupied and all of the 2NN sites 
are occupied. For R\ = R2 = —1/2, the superantiferro- 
magnetic configuration has less NN and 2NN repulsion, 
thus it has lower free energy. We can expect that if the 
NN repulsion is sufficiently stronger (Ri < 2i? 2 ), the sta- 
ble phase in the middle would be antiferromagnetic. 



V. CONCLUSION 

In the same spirit that the nearest neighbor lattice gas 
helps to understand the liquid-gas transition, we have 
used lattice models with an extended bonding range to 
understand how the formation of interpenetrating open 
networks can give rise to liquid- liquid transitions. The 
specific mechanism of interpenetration seems most appli- 
cable to recently studied DNA functionalized nanoparti- 
cles, where the range of bonding can be quite large com- 
pared with the core exclusion |27j . 

The presence of lattice sites in the lattice models we 
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have studied provides a predetermined regular frame- 
work for the interpenetration. As a result, it is natu- 
ral to generate distinct phases on the sublattices defined 
by second neighbors (or more complex sublattices in the 
case of the third neighbor model). However, in contin- 
uum systems, there is no such underlying lattice struc- 
ture. For systems with highly directional interactions, 
a network structure may appear naturally. Interpene- 
tration for tetrahedral networks is facilitated by the fact 
that the empty spaces of the tetrahedral network make 
a complementary tetrahedral network; the same is true 
for cubic network. For such systems, the free energy of 
distinct networks can be significantly lower than that of 
a distorted, connected structure. As a result, the system 
will preferentially phase segregate at densities where dis- 
tinct networks with few defects are not possible. Thus, 
for networked structures lacking such symmetry, the in- 
terpenetration could be frustrated, thereby eliminating 
the free energy gap between distorted networks and dis- 
tinct interpenetrating networks. Without such a free en- 
ergy gap, no phase separation will occur. In particular, 
spherically symmetric systems with a long-ranged bond- 
ing term would likely not result in multiple phases due 
to interpenetration. However, packing constraints may 
still give rise to polyamorphic behavior in step potentials 
with carefully chosen step sizes [2U1 [HJ [22] . 

The lattice models we have considered might be made 
more specific to the problem of DNA functionalized 
nanoparticles by including specific interactions that could 
mimic the molecular recognition of DNA, and by control- 
ling the number of bonds that a given site can participate 
in. Such a model might be useful for developing a qual- 
itative understanding of how mixing different species of 
DNA functionalized nanoparticles might give rise to net- 
works which are chemically distinct, and in their thermo- 
dynamic properties. Additionally, constraining the num- 
ber of neighbors in simple spherical potentials is known 
to dramatically alter the phase behavior [51, 52 551 154] . 
Thus it can be expected that changing the number of 
DNA strands attached to a core nanoparticle could have 
a similar dramatic effect. 
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APPENDIX 

During our initial study of lattice gas models with sec- 
ond neighbor interactions, we attempted to reproduce 
results of Binder and Landau for the second neighbor 
antiferromagnetic case £2/^1 = 1/2 and e\ < 0. Our re- 
sults confirmed the major finding of the earlier work |32j . 
Due to the dramatic increase in computing power over 28 
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FIG. 9: Mean-field phase diagram for the case £2/61 = 1/2, 
ei > 0, as a correction for Fig. 6(b) in Ref. [32] . All transi- 
tions are first order. 



years since Ref. (32] was published, we could examine the 
phase behavior much more finely. As a result, we found 
that Binder and Landau identified second order transi- 
tions which are in fact actually first order. Specifically, 
the highest and lowest density transitions were misiden- 
tified in Ref. [32] . We provide a corrected version of this 
phase diagram in Fig. [9] (compared with Fig. 6(b) in 
Ref. (35]). The jump in the slope of free energy is too 
small to be detected with the accuracy of computation 
at that time. 
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